function create_data(filename, job, options)
    if strcmpi(job, 'margin')
% create_data('data/margin2', 'margin', struct('N',100,'dimension',2,'min_dist',1,'max_count',100,'sparse',1))
        data = margin_data(options);
        save(filename, 'data');
    elseif strcmpi(job, 'plot')
        load(filename);
        plot_data(data);
    elseif strcmpi(job, 'plot_decision')
% create_data('data/margin2', 'plot_decision', model)        
        load(filename);
        plot_decision(data, options);
    elseif strcmpi(job, 'plot_regpath')
% create_data([], 'plot_regpath', struct('model',model,'w_name','w','C_name','C'))
        plot_regpath(options);
    end
end

function data = margin_data(options)
    N = options.N;
    d = options.dimension;
    min_d = options.min_dist;
    max_count = options.max_count;
    sparse = options.sparse;
    
    w = randn(1, d)*5; 
    if ~strcmpi(sparse, 'default')
        w(sparse+1:end) = 0;
    end
    b = randn;
    norm_w = norm(w,2);
    
    X = zeros(N, d);
    y = zeros(N, 1);
    for i = 1:N
        count = 0;
        while true
            x = randn(1, d)*5;
            f = dot(w, x) + b;
            dist = abs(f) / norm_w;
            if (dist >=min_d) 
                X(i,:) = x;
                y(i) = sign(f);
                break;
            end
            count = count + 1;
            if (count > max_count)
                error('Cannot make new points');
            end
        end
    end
    
    data.w = w;
    data.b = b;
    data.X = X;
    data.y = y;
end

function plot_data(data)
    red = data.y == 1;
    plot(data.X(red,1), data.X(red, 2), 'r+');
    hold on
    plot(data.X(~red,1), data.X(~red, 2), 'bo');
    hold off
end

function h = plot_decision(data, options)
    xrange = [min(data.X(:,1))-1 max(data.X(:,1))+1];
    yrange = [min(data.X(:,2))-1 max(data.X(:,2))+1];
    h = show_boundary(xrange, yrange, @(XX)linear_boundary(XX, options.w, options.b));
    
    red = data.y == 1;
    plot(data.X(red,1), data.X(red, 2), 'r+');
    plot(data.X(~red,1), data.X(~red, 2), 'bo');
    
    legend('Class +1', 'Class -1','Location','NorthOutside', ...
            'Orientation', 'horizontal');
    xlabel('x');
    ylabel('y');
end

function idx = linear_boundary(XX, w, b)
    idx = (XX * w + b > 0) + 1;
end

function h = show_boundary(xrange, yrange, boundary_function)
    inc = 0.1;
    [x, y] = meshgrid(xrange(1):inc:xrange(2), yrange(1):inc:yrange(2));
    image_size = size(x);
    XX = [x(:) y(:)];
    idx = boundary_function(XX);
    decisionmap = reshape(idx, image_size);
    
    h = clf;
    imagesc(xrange,yrange,decisionmap);
    hold on;
    set(gca,'ydir','normal');
    cmap = [1 0.8 0.8; 0.95 1 0.95; 0.9 0.9 1];
    colormap(cmap);
    
end

function plot_regpath(options)
    model = options.model;
    w_name = options.w_name;
    C_name = options.C_name;
    
    N = length(model);
    d = eval(['length(model(1).' w_name ')']);
    
    w = zeros(d,N);
    C = zeros(1,N);
    
    for i = 1:N
        for j = 1:d
            w(j,i) = eval(['model(i).' w_name '(j)']);
        end
        C(i) = eval(['model(i).' C_name]);
    end
    
    h = clf;
    semilogx(C, w);
end